\chapter{万有引力定律推导 Law of Universal Gravitation}
\section{万有引力定律}
\subsection{导出万有引力定律}
现在我们要利用前面的推导结论，推导万有引力定律。

根据式 \ref{kepler3rdlaw26} 和式 \ref{Elliptic_equation10}，可得

\begin{equation}
	\label{kepler3rdlaw29}
	\frac{C^2}{p}=\frac{4\pi^2a^3}{T^2}=k^2
\end{equation}

联立 \ref{bineteqnewton11_1}和 \ref{Elliptic_equation_polar}，解得

\begin{equation}
	\label{kepler3rdlaw30}
	u=\frac{1+ecos\theta}{p}
\end{equation}

上式对角度求导，得

\begin{equation}
	\label{kepler3rdlaw31}
	\frac{du}{d\theta}=\frac{-esin\theta}{p}
\end{equation}

再对角度求导，得

\begin{equation}
	\label{kepler3rdlaw32}
	\frac{d^2u}{d\theta^2}=\frac{-ecos\theta}{p}
\end{equation}

联立式 \ref{kepler3rdlaw32}和 \ref{kepler3rdlaw30} 、 \ref{bineteqstd}，解得

\begin{equation}
	\label{kepler3rdlaw33}
	-\frac{F}{m}=C^2u^2(\frac{d^2u}{d\theta^2}+u)=C^2u^2(\frac{-ecos\theta}{p}+\frac{1+ecos\theta}{p})=\frac{C^2}{p}u^2=k^2\frac{1}{r^2}
\end{equation}

即

\begin{equation}
	\label{kepler3rdlaw34}
	F=-\frac{mk^2}{r^2}
\end{equation}

根据牛顿第三定律，A吸引B的力与B吸引A的力是一对作用力和反作用力，它们应当有相同的性质，即B对A的吸引力的表达式为

\begin{equation}
	\label{kepler3rdlaw35}
	F'=-\frac{Mk'^2}{r^2}
\end{equation}

由牛顿第三定律知道，F与F'大小相等，即

$F=F'$

联立上面三式，得

\begin{equation}
	\label{kepler3rdlaw36}
	-\frac{mk^2}{r^2}=-\frac{Mk'^2}{r^2}
\end{equation}

移项后得

\begin{equation}
	\label{kepler3rdlaw37}
	\frac{k^2}{M}=-\frac{k'^2}{m'}
\end{equation}

因为k和M是已知量，所以可以令

\begin{equation}
	\label{kepler3rdlaw38}
	G=\frac{k^2}{M}
\end{equation}

即
\begin{equation}
	\label{kepler3rdlaw45}
	k^2=GM
\end{equation}

上式代入式 \ref{kepler3rdlaw34}得

\begin{equation}
	\label{kepler3rdlaw46}
	F=-\frac{GMm}{r^2}
\end{equation}

这就是万有引力定律。

上面是根据椭圆轨道方程推导得到万有引力定律的。如果是双曲线轨道或抛物线轨道，是否也能导出万有引力定律呢？

因为圆锥曲线极坐标形式一般方程为 \ref{conic_equation_polar}，它比椭圆极坐标方程 \ref{Elliptic_equation_polar}只在右端分子多了系数e，从而u及其1、2阶对角度的导数在形式上只比 \ref{kepler3rdlaw30}、 \ref{kepler3rdlaw31}、 \ref{kepler3rdlaw32}分母多一个系数e，因此在式  \ref{kepler3rdlaw33}第2个等号后也只是分母多一个系数e，因此不改变  \ref{kepler3rdlaw34}的形式。这就是说，只要满足圆锥曲线的极坐标方程，就可以得到万有引力定律。这个结论是否合理呢？

这个结论是合理的。万有引力可以使得质点运动到平面内任何位置。在二体运动中，质点的轨道是圆锥曲线，根据不同的初速度，轨道曲线可以是点、直线、圆、椭圆、抛物线、双曲线。质点运动只有在圆或椭圆轨道是稳定的，在其它类型轨道上运动将只有经过焦点附近时才能被短暂地观测到。稳定的椭圆轨道是质点的基态，质点在得到外部能量以后，会转移到能量更高的椭圆轨道运行。如果质点能量达到临界能量，质点会进入抛物线轨道。如果质点能量超过临界能量，质点将进入双曲线轨道，永远弹射出去离开系统，不再回归。这种现象称为粒子碰撞，粒子对撞机就是根据类似原理设计的。

\begin{equation}
	\label{kepler3rdlaw47}
	F=-\frac{GMm}{r^2}
\end{equation}

\begin{equation}
	\label{Elliptic_equation_polar}
	r=\frac{p}{1+ecos\theta}
\end{equation}

\subsection{导出行星轨道定律}
现在我们推导行星轨道定律。

行星绕太阳运行，受到万有引力：

\begin{equation}
	\label{planetorbit1}
	F=-\frac{GMm}{r^2}=-GMmu^2
\end{equation}

代入比耐方程 \ref{bineteqstd} 得到

\begin{equation}
	\label{planetorbit2}
	h^2u^2(\frac{d^2u}{d\theta^2}+u)=GMu^2
\end{equation}

化简为

\begin{equation}
	\label{planetorbit3}
	\frac{d^2u}{d\theta^2}+u=\frac{GM}{h^2}
\end{equation}
令
\begin{align}
	v&=u-\frac{GM}{h^2}\label{planetorbit31}\\
	\frac{d^2v}{d\theta^2}&+v=0\label{planetorbit32}
\end{align}

式\ref{planetorbit32}是简谐振动方程，通解为三角函数，解为

\begin{align}
	v&=Acos(\omega \theta-\theta_0)\label{planetorbit41}\\
	\omega^2&=\frac{k}{m}=1\label{planetorbit42}\\
	u&=v+\frac{GM}{h^2}\\
	u&=\frac{GM}{h^2}+Acos( \theta-\theta_0)\label{planetorbit4}
\end{align}

代入r=1/u，得

\begin{align}
	r&=\frac{1}{\frac{GM}{h^2}+Acos( \theta-\theta_0)}\label{planetorbit5}\\
	r&=\frac{\frac{h^2}{GM}}{1+\frac{Ah^2}{GM}cos( \theta-\theta_0)}\label{planetorbit6}
\end{align}
振幅A 和初相位$\theta_0$为积分常数，由初始条件定：$t = 0,  \theta=\theta_0$.

上式即为轨道方程，几何上为极坐标表示的圆锥曲线方程。

对比圆锥曲线的极坐标方程
\begin{equation}
	\label{Elliptic_equation_polar}
	r=\frac{p}{1+ecos\theta}
\end{equation}
可以看出
\begin{align}
	\mbox{偏心率}e&=\frac{Ah^2}{GM}\label{planetorbit7}\\
	\mbox{半正焦弦长}p&=\frac{h^2}{GM}\label{planetorbit8}
\end{align}

因此，行星在万有引力作用下绕太阳运行的轨道是圆锥曲线，稳定的轨道是椭圆曲线，这就是开普勒发现的行星运动第一定律。

至此，我们根据比耐方程，从开普勒行星运动定律推导出了万有引力定律，也从万有引力定律推导出了开普勒行星运动定律。可以证明这两个定律反映了物质世界的普遍规律，从宏观天体到基本粒子都适用这两个定律。

后面我们将要使用它们推导库仑定律和电磁场方程以及其它物理定律。
\subsection{偏心率e与能量E对与圆锥曲线轨道的影响}
圆锥曲线\ref{Elliptic_equation_polar}上某点的轨迹与偏心率e有如下关系：

1) e=0，轨迹为一点或一个圆；

2) 0<e<1，轨迹为椭圆；

3) e=1(即到P与到L距离相同)，轨迹为抛物线；

4) e>1，轨迹为双曲线。

知道了偏心率e就可以判定圆锥曲线的形状，而e需根据初始条件来确定，e是轨道的几何常数，在力学中，希望借助动力常数如能量等来作为轨道类别的判据。

由于距离平方反比引力的形式为：
\begin{align}
	\vec{F(r)}&=-\frac{k^2}{r^2}\vec{e_r}
\end{align}
若取无穷远处的势能为零，则质点在距力心r处的势能为：
\begin{align}
	V(r)&=\int-F(r)dr=\int_\infty^r\frac{k^2}{r^2}dr=-\frac{k^2}{r}\\
	\mbox{又}\dot{r}&=\frac{dr}{dt}=\frac{dr}{d\theta}\frac{d\theta}{dt}=\dot{\theta}\frac{dr}{d\theta}=(r^2\dot{\theta})\frac{1}{r^2}\frac{dr}{d\theta}=\frac{h}{r^2}\frac{dr}{d\theta}\\
	\mbox{机械能守恒：}T+V(r)&=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)-\frac{k^2}{r}=E\label{mechanicalenergyconservation}\\
	T+V(r)&=\frac{1}{2}m[\frac{h^2}{r^4}(\frac{dr}{d\theta})^2+\frac{h^2}{r^2}-\frac{2k^2}{mr}]=E\\
\end{align}
解出$\frac{dr}{d\theta}$，并分离变量，可得：
\begin{align}
	\frac{h^2}{r^4}(\frac{dr}{d\theta})^2&+\frac{h^2}{r^2}-\frac{2k^2}{mr}-\frac{2E}{m}=0\\
	(\frac{h}{r^2}\frac{dr}{d\theta})^2&=-\frac{h^2}{r^2}+\frac{2k^2}{mr}+\frac{2E}{m}\\
	\frac{hdr}{r\sqrt{\frac{2E}{m}r^2+\frac{2k^2r}{m}-h^2}}&=d\theta\\
	\mbox{利用积分公式：}\int\frac{dx}{x\sqrt{a+bx+cx^2}}&=-\frac{1}{\sqrt{-a}}\arcsin\frac{bx+2a}{x\sqrt{b^2-4ac}},(a<0)\\
	\mbox{两边分别积分，可得：}\arcsin\frac{2k^2r-2h^2m}{r\sqrt{4k^4+8Eh^2m}}&=\theta+\frac{3\pi}{2}-\theta_0\\
	\mbox{故}\frac{2k^2r-2h^2m}{r\sqrt{4k^4+8Eh^2m}}&=\sin(\theta+\frac{3\pi}{2}-\theta_0)\\
	2k^2-\frac{2h^2m}{r}&=-\sqrt{4k^4+8Eh^2m}\cos(\theta-\theta_0)\\
	\frac{2h^2m}{r}&=2k^2+\sqrt{4k^4+8Eh^2m}\cos(\theta-\theta_0)\\
	r&=\frac{2h^2m}{2k^2+\sqrt{4k^4+8Eh^2m}\cos(\theta-\theta_0)}\\
	r&=\frac{h^2m/k^2}{1+\sqrt{1+2Eh^2m/k^4}\cos(\theta-\theta_0)}\label{conicsectionenergyEq}\\
	\mbox{比较圆锥曲线方程}r&=\frac{p}{1+e\cos(\theta-\theta_0)}\\
	\mbox{知}e&=\sqrt{1+\frac{2mh^2}{k^4}E}\\
	\mbox{两边平方，得 }e^2&=1+\frac{2mh^2}{k^4}E\\
	\mbox{因为k为高斯引力常数，满足 }k^2&=GM\\
	\mbox{所以 }e^2&=1+\frac{2mh^2}{(GM)^2}E\\
\end{align}
这个结果也可用机械能守恒和角动量守恒来求，但物理意义远不及这个积分丰富，因为积分中直接含E。

机械能守恒有方程\ref{mechanicalenergyconservation}

由轨道方程$r=\frac{p}{1+e\cos\theta}$知，当$\theta=0$时，质点与力心距离最近，用$r_{min}$表示，且
\begin{align}
	r_{min}&=\frac{p}{1+e}=\frac{mh^2}{k^2(1+e)}
\end{align}
由极值条件知$\dot{r}=0$,即此处径向速度为0。

在$r_{min}$处的总能量可以表示为：
\begin{align}
	E&=\frac{mh^2}{2r^2_{min}}-\frac{k^2}{r_{min}}=\frac{1}{2r^2_{min}}(mh^2-2k^2r_{min})=(\frac{k^2(1+e)}{2mh^2})^2(mh^2-2k^2\frac{mh^2}{k^2(1+e)})\\
	&=(\frac{k^2(1+e)}{2mh^2})^2(mh^2)(1-\frac{2k^2}{k^2(1+e)})=(\frac{k^2(1+e)}{2mh^2}\frac{k^2(1+e)}{k^2(1+e)})(k^2(1+e)-2k^2)\\
	&=\frac{k^4}{2mh^2}(e^2-1)
\end{align}
由总能量
\begin{align}
	E&=\frac{k^4}{2mh^2}(e^2-1)\label{RelationEe1}
\end{align}
反解e跟E的关系为：
\begin{align}
	e&=\sqrt{1+\frac{2mh^2}{k^4}E}\label{RelationEe2}\\
	\mbox{或}e&=\sqrt{1+\frac{2mL_0^2}{mk^4}E}\label{RelationEe3}
\end{align}


由e和E的关系可以看出：

0、E<0(e=0)-------圆

1、E<0(e<1)-------椭圆

2、E=0(e=1)-------抛物线

3、E>0(e>1)-------双曲线

在椭圆中，由近日点$r_{min}(\theta =0) \mbox{和远日点}r_{max} (\theta  = \pi  )$ 

以及椭圆的定义知：
\begin{align}
	2a&=r_{min}+r_{max}=\frac{2p}{1-e^2}=\frac{-k^2}{E}\label{RelationEe4}
\end{align}
\subsubsection{椭圆轨道情况}
如果球或粒子运行在椭圆形轨道，根据式\ref{planetorbit7}，必须0<e<1，因此必须
\begin{align}
	0&<\frac{Ah^2}{GM}<1\label{planetorbit81}\\
	0&<A<\frac{GM}{h^2}\label{planetorbit82}\\
\end{align}
椭圆轨道满足振幅A>0，$0<\frac{Ah^2}{GM}<1$。

\subsubsection{圆轨道情况}
如果球或粒子运行在圆形轨道，根据式\ref{planetorbit7}，必须e=0，因此必须$M\rightarrow\infty$或A=0。

振幅A为二体系统质心的运动轨迹曲线的半径。

A=0表示振动方程的振幅为0，振幅为0表示球或粒子绕质心公转的线速率变化为0，即线速率变化幅度已经缩小为一个在圆心的点。

\subsubsection{轨道的稳定性}
研究轨道的稳定性对航空航天领域具有十分重要的意义。

我们讨论质点在受到有心力作用时的圆形轨道问题。

设有心力势场是V(r)，则有心力为：
\begin{align}
	\vec{F}=F(r)\vec{e_r}=-\nabla V(r)&=-\frac{dV}{dr}\vec{e_r}\label{circularorbitstability1}\\
	\mbox{有心运动中，机械能守恒：}E&=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)+V(r)=\frac{1}{2}m\dot{r}^2+\frac{mh}{2r^2}+V(r)=\frac{1}{2}m\dot{r}^2+U(r) \label{circularorbitstability2}\\
	U(r)&=\frac{mh}{2r^2}+V(r) \label{circularorbitstability3}
\end{align}
它等效于质点沿径向作一维运动时的机械能守恒。

我们知道，要使质点在势场中保持稳定，只有处于势场中势能最小处，由数学知识可知，只需 U(r) 的一阶导数等于0，二阶导数大于0，即：
\begin{align}
	\frac{dU(r)}{dr}&=0\\
	\frac{d^2U(r)}{dr^2}&>0\\
	\frac{dU(r)}{dr}&=-\frac{mh^2}{r^3}+\frac{dV(r)}{dr}=-\frac{mh^2}{r^3}-F(r)=0\label{circularorbitstability4}\\
	\frac{d^2U(r)}{dr^2}&=\frac{3mh^2}{r^4}-F'(r)>0\label{circularorbitstability5}
\end{align}
若解出此时的$r=r_0$，显然稳定条件为：
\begin{align}
	F(r_0)&=-\frac{mh^2}{r_0^3}\label{circularorbitstability6}\\
	F'(r_0)&>\frac{3mh^2}{r_0^4}\label{circularorbitstability7}\\
	\frac{F'(r_0)}{F(r_0)}&>\frac{-3}{r_0}\label{circularorbitstability8}
\end{align}
为具体讨论，不妨设引力的大小与距离成n次幂关系：
\begin{align}
	F(r)&=-\frac{k^2}{r^n}\label{circularorbitstability9}\\
	F'(r)&=\frac{nk^2}{r^{n+1}}\label{circularorbitstability10}\\
	\frac{F'(r)}{F(r)}&=\frac{\frac{nk^2}{r^{n+1}}}{-\frac{k^2}{r^n}}=-\frac{n}{r} >\frac{-3}{r}\label{circularorbitstability11}\\
	n&<3
\end{align}
可见圆轨道的稳定性的条件是n<3。

这说明质点在受到距离平方反比引力作用时，圆轨道是稳定的。

注：式\ref{circularorbitstability9}、\ref{circularorbitstability10}、\ref{circularorbitstability11}假定和推导在数学上是不严格的。如果假定了式\ref{circularorbitstability10}，则必须重新推导能量公式及式\ref{circularorbitstability5}。n<3是圆和椭圆轨道稳定的必要条件。n=2是圆和椭圆轨道稳定的充分必要条件。
